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We study the coalescence of two drops of an ideal fluid driven by surface tension. The 
velocity of approach is taken to be zero and the dynamical effect of the outer fluid 
(usually air) is neglected. Our approximation is expected to be valid on scales larger than 
l v — pv 2 /a, which is lOnm for water. Using a high-precision boundary integral method, 
we show that the walls of the thin retracting sheet of air between the drops reconnect in 
finite time to form a toroidal enclosure. After the initial reconnection, retraction starts 
again, leading to a rapid sequence of enclosures. Averaging over the discrete events, we 
find the minimum radius of the liquid bridge connecting the two drops to scale like 
r b cx t 1 ' 2 . 



1. Introduction 



Drop coalescence arises in many different contexts, and is crucial to our understanding 
of free surface flows ( Eggers| |l997| ). Examples are printing applications ( Chaudhary fe 



Maxworthyl |1980| ; |Wallace| |200l| ), drop impact on a fluid surface (pguz fc Prosperetti| 



1990|), and the coarsening of drop clouds and dispersions ( 


MacPhee et al. 


2002 


; Jury 



2 

et al. 



L. Duchemin, J. Eggers and C. Josserand 



1999 



Verdiei 200C). After the two surfaces have merged on a microscopic scale, 



surface tension drives an extremely rapid motion, usually impossible to resolve in either 



experiment (Bradley & Stow 1978; Mcnchaca-Rocha et al. 2001) or simulation (Lafaurie 



et al. 1994). Thus theory is needed to investigate a possible dependence on initial condi- 



tions, development of small-scale structures during merging, and to estimate the typical 
time required for merging. 

A large body of work exists on this problem in the case that viscosity is dominant 
and the motion is described by Stokes' equation. In the absence of an outer phase this is 



known as the "viscous sintering problem" (Frcnkcl 1945; Hopper 1993; Martinez-Herrera 



& Derby 1995), the inclusion of an outer phase is important for many problems governing 



the coarsening of dispersion (Nikolayev et al. 1996; Verdiei 2000). For the two-dimensional 



problem (i.e. for the merging of cylinders) exact solutions exist (Hopper 1990; Richardson 
1992| ; prowdy| |2002| , [To appear| ), which were shown (Eggers et al\ |1999|) to be asymp- 



totically equivalent to their three-dimensional counterparts. The presence of an outer 



fluid leads to the formation of a toroidal bubble during merging (Eggers et al. 1999), 
significantly modifying the dynamics. 

Fig. |l| shows two equal drops of radius R being connected by a liquid bridge of radius 
rt>, which is rapidly being pulled up by surface tension. The local Reynolds number of 
this flow can be estimated as Re = ari } /(pv 2 ) J where a is the surface tension, p the 
density, and v the kinematic viscosity. Thus, regardless of the value of the viscosity, the 
Reynolds number is always small in the initial phases of the merging, which is equivalent 
to demanding that rj, <C t v , where l v = v 2 p/a is the viscous length scale. However, t v is 



often very small (140 A for water, and 4 A for mercury (Eggen 1997)), so r?, ^ i v for 
a large part of the evolution, and inviscid theory can be applied. Thus for a wide range 
of practical problems the almost inviscid regime, which is the topic of this letter, is the 
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Figure 1. Initial condition. Two drops touching at a point are joined by a liquid bridge of 
radius »v The inset shows the width of the gap just above the meniscus to be w = r^. The gap's 
walls are nearly straight on the scale of w. 

most relevant. Typically, the viscous regime will serve as an inner layer that defines the 
initial condition for the inviscid problem we are interested in. In general, we do not have 



to worry about the initial process of reconnection (Amarouchenc et al. 2001), which for 
clean fluids is expected to take place over a microscopically small area. 

In the case of a head-on collision of two drops with relative velocity V, considered 



in ( Oguz fc Prosperetti 1989 ), a purely geometrical consideration predicts ss \fVRt 
for two overlapping circles. The corresponding speed of merging is of the same order 
as the surface-tension-driven merging to be described below, so V has thus to be taken 
into account. However, we will restrict ourselves here to the case where the velocity 
of approach is vanishingly small, a condition that is easily realizable experimentally 
( Menchaca-Rocha et al. [2001 ) . We also do not treat the dynamical effect of an outer fluid 
like air, which might become important as the lubrication layer between the approaching 



4 L. Duchemin, J. Eggers and C. Josserand 

drops becomes very thin (Eggers et al. 1999; Yiantsos & Davie 1991). However, this 
approximation is consistent with the assumption of a small velocity of approach. 



2. Initial conditions and scaling laws 

We consider two identical drops of radius R touching at a point where a thin liquid 
bridge of size r& connects the two drops initially (cf. figure |l|). The general problem of 
drops of different radii only changes a prefactor in the gap width between the drops 
( Eggers] 199S). For the inviscid dynamics considered here, all parameters of the problem 



can be scaled out by writing the time and space coordinates in units of yj pR 3 ja and R, 
respectively. Assuming that the vorticity generated by the initial viscous motion can be 
neglected, and using incompressibility, the velocity potential (p obeys 

Aip = 0. (2.1) 

The boundary condition on the free surface amounts to a balance between surface tension 
and Bernoulli pressures flOguz fc Prosperctti||l989| ): 



|f + ^(V^) 2 - K = 0, (2.2) 



where k is the mean curvature of the interface. 



We have to solve (2T), (2/2) with the initial condition shown in figure [l| assuming that 
the bridge radius r b is initially very small (typically 10 in our numerical simulations). 
Away from the point of contact at z = 0, but for h <C 1 the surface has the form 
h(z) = (2Z) 1 / 2 and h(z) = (~2z) 1 ^ 2 for z > and z < 0, respectively. The width of the 
gap at a height r is thus 

w = r 2 (r b < r < 1) (2.3) 
and since dw/dr <C 1, the walls are nearly parallel. Thus the meniscus, which owing to 
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radial symmetry is located along a ring of radius n, is being pulled straight up by a force 

2a per unit length. 



Assuming that the profile in region (2.3) matches onto the bridge on the scale r ~ r&, 



the curvature at the meniscus can be estimated as Kb w r, 2 , much larger than the axial 



curvature r^ 1 of the liquid bridge. Thus, as already argued in ( Eggcrs et al. 1999 ), the 
axial curvature can be neglected for rj < 1 and the problem becomes effectively two- 
dimensional, equivalent to the merging of two fluid cylinders. Thus a model problem 



( |Oguz fc Prospcrctti 1989 ; Eggers 1998| ) for the initial motion of the meniscus is that of 



the two dimensional, straight slot shown in the inset of figure [l]. The eventual widening 
of the gap can be neglected on the scale of the gap width w. 

The results of our computations for the full three-dimensional problem, to be explained 
in more detail below, are shown in figure ^. As the meniscus retracts, the rapid fluid flow 
past the sides of the gap creates an under-pressure as described by Bernoulli's equation 
( |2.2| ), which in turn causes the end to expand into a bubble. As the bubble increases 
in size, capillary waves are excited in its wake, with amplitude roughly proportional to 
the bubble radius. Thus after the amplitude of the capillary wave has grown to the half 
width of the slot w/2, its two sides touch and reconnect at a time r c . Since the width is 
the only length scale in the problem, it follows that the total length r c the meniscus has 
retracted up to the point of reconnection is proportional to w, while the time r c required 
scales like w 3 ^ 2 . We thus have 

r c = r w, t c = t w 3/2 , (2.4) 

where ro, To are constants to be determined numerically. Below we find in fact ro = 
10, r = 7.6. 

After the two sides of the gap have reconnected, this new initial condition looks very 
similar to the original one, except for a non-trivial velocity field that remains. But since 
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Figure 2. A sequence of profiles showing the retraction of the initial meniscus for rb(0) = 1CP 5 . 
At a time r c = to w 3 ^ 2 = 7.6 w 3 ^ 2 the walls of the gap touch and the minimum radius z m i„ goes 
to zero. The distance of this point from the initial tip of the meniscus is r c = row = 10 w. 

most of the resistance to the motion before reconnection is due to the large bubble that 
was left behind, this velocity can be neglected relative to the velocity to be generated at 
the next stage of the motion (more detailed estimates are given below) . This means that 
at each step the same motion repeats itself, but with a slightly larger radius r^. At the 



n-th step we can thus write, analogous to ( Eggers |199q ) 



r n+1 — r n — r — r„(r n \ 2 

r b 7 b — 7 c — r (r b ) , 



and for the times t n of successive pinching events: 



tn+l - t» = T (r^) 3 . 

For very small initial r\> reconnection occurs in rapid succession, with small relative 
change of the variables. We can thus write n, as a smooth function of t, obeying the 



Inviscid coalescence of drops 



7 



differential equation 



dr b 
dt 



r 1 



T n 



(2.5) 



which gives, after integration : 




(2.6) 



The scaling law ( |2.6| ) is the central result of the present letter. Eventually, when is of 
the same order than the drop radius, the widening of the channel overcomes the growth 
of capillary waves, and the enclosure of bubbles stops. This is when the time scale of 
retraction r ss r^ro/(2ro) is shorter than r c w r^rf characterizing reconnection. Thus 
reconnection will cease when > l/(2ro) = 0.05. We have determined numerically that 
no more voids are entrapped for r& > 0.035, in good agreement with our theoretical 
estimate. Below we present detailed numerical tests of the scaling predictions, and inves- 
tigate further the crucial stage of bubble growth, from which we are able to extract the 
numerical constants ro, tq. 

3. Boundary integral method 

If the flow can be considered potential and incompressible, the use of a boundary 
integral method is advantageous, since the velocity field can be calculated from the 
interface shape. Thus one only needs to keep track of the interface, represented by a 
one-dimensional curve, and grid refinement can be done very efficiently. The majority 
of these boundary integral methods require smoothing of the surface, in order to avoid 
short wave length instabilities. The method briefly presented here does not require any 
explicit smoothing, except for a redistribution of the points around the tip at every time 
step. This redistribution can act as a smoothing, but no damping of instabilities, such 
as an artificial surface viscosity, has been used. 

The dipole formulation used here is very close to the one described by Baker, Meiron 



L. Duchemin, J. Eggers and C. Josserand 



and Orszag (Baker et al. 1980), but it needs to be refined to be able to resolve the very 
disparate scales of the drops and of the highly curved region close to the meniscus. At a 
given time step, we expect the velocity potential tp to be known, from which we calculate 
the normal and the tangential velocity of the surface. This velocity is then used to advect 



the surface, and to advance (p using Bernoulli's equation (2.2). The tangential velocity is 
calculated directly by differentiating with respect to the arclength along the interface: 



Ui 



d(p 
1h 



(3.1) 



to compute the normal component, we use the vector potential A of the velocity field, 
u = V x A: 

1 drAg 



r ds 



Following (Baker et al. 1980), we first compute the dipole density (i from 



<p(M) = fx(M) 



-in 



(M(M) - MM')) 



d (i 



dS', 



(3.2) 



(3.3) 



dn \X / 

where A is the distance between points M and M' on the surface. The appearance of 
(J>(M) in the integrand serves to subtract the singularity of the normal derivative. Once 
ji is known, it can be used to calculate the vector potential: 



MM) = ^f (M(M') - n{M)) n x V s (j)dS'. 



(3.4) 



Classical iterative solutions of (^3), (3^) were found to fail for very small bridge radii, 
so ( |3.3| ),( |3T| ) were solved by matrix inversion instead. A simple trapezoidal rule was used 
to convert the equations into linear systems, which was then solved by LU decomposi- 
tion. In order to compute the curvature of the surface and the tangential derivatives in 



( p.l[ ),(3.2), we re-parametrized the integrals by introducing a new integration variable 
which equals i at grid-point i. This avoids instabilities in the cubic spline interpolation 
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that would otherwise be present if two points come very close together, as it happens at 

the tip. 

At each time step, the Bernoulli equation and the kinematic condition were used to 



advance the solution using a Crank-Nicolson scheme (Press et at 1992). The implicit 
equations were solved by iteration, which required less than 10 iterations until a relative 
error of 10 -5 in the velocity potential was reached. An explicit Runge-Kutta fourth order 
scheme was also tested, but found to be too unstable for small values of r&. 

We also redistribute grid-points at every time step according to the their distance from 
the tip. Cubic splines are used to interpolate to the new points. At each time step points 
are placed on the free surface with grid spacing 5; typical values are shown in figure ||. 
This spacing is used up to a distance of 40 r\ from the tip, after which it is gradually 
increased in steps of 2, since much lower resolution is required far from the tip. 

4. Reconnection 

As we have explained above, the retraction of the meniscus is interrupted by the 
reconnection of the two sides of the gap, and the distance r c by which the meniscus 
recoils as well as the time r c required is given by the scaling relations (2.4). In figure || we 



define typical quantities characterizing the retraction of the meniscus. The minimum gap 
radius z m i n marks the first trough of a train of capillary waves that is generated by the 



growing bubble. Note that in the corresponding simulation in (Oguz & Prosperetti 1990) 
(cf. figure 4) there is little or no indication of this growth of capillary waves. We suspect 
that these authors did not follow the retraction for sufficiently long times, and that the 
low resolution of their simulation introduced additional damping, which smoothed out 
the capillary waves. 

As seen in figure the time dependence of the minimum gap radius z m i n converges 
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Figure 3. The minimum gap radius z m i n /r%(Q) plotted against t/rf (0). The initial value of rj, 
is 10 -5 and the three resolutions correspond to the minimum distance between points in the tip 
region. The linear extrapolation gives tj,j no ft — 7.6 r%(0) 

towards a close to linear behavior as the resolution is increased. Extrapolation towards 
z-min — thus gives a reliable estimate of the time required for reconnection. Although 
the walls of the gap do not interact physically, errors of our boundary integral description 
grow large as two surfaces become close to each other. The reason is that the distance A 
between points varies on scale z m j„ close to the minimum, so the grid spacing 6 always 
needs to be smaller than z m i n . 

From the simulations we deduce the values = 10 and tq = 7.6 for the reduced 
retraction length and time already reported in section 2. Here the underlying assumption 
is that the dynamics is controlled by the local gap width alone. To test this idea, we have 
computed a sequence of pinch events as shown in figure ||. When z m i n has gone down to 
about 10 % of the local gap radius w/2, the gap is cut at about w/2 behind the minimum 
and new points are introduced along the new surface. Our method of redistributing points 
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Figure 4. Successive entrapment of voids during the coalescence for an initial liquid bridge 
radius of r& = 0.008. After every reconnection, the void is extracted from the profile and a new 
computation begins, with a null initial velocity field. 

automatically introduced a certain smoothing, which was enough for the simulation to 
continue. Obtaining a new initial condition for the velocity profile proved to be much 
more difficult. Simply extrapolating the velocity potential tp before the surgery to the 
new initial condition led to instabilities that could no longer be controlled numerically, so 
instead we had to put the velocity field to zero. This is justified by the fact that the gap 
position very quickly re-assumes its retraction velocity after the bubble is left behind, as 
we discuss in more detail below. 

As illustrated in figure [|, this leads to a self-similar succession of pinch-off events. Each 
simulation was started from a new value of the bridge radius r£ . The typical gap width 
at the meniscus is then w — (r£) 2 . A more quantitative test of the scalings employed 
in section 2 is presented in figure |^, where we plot the bridge radius rf, as a function of 
time and, in the inset, r c /r c = (to/tq)/^ as function of the bridge radius at the time 
of pinching. The excellent agreement with the predicted scaling behavior confirms our 
assumption that the local dynamics only depends on the gap width at the corresponding 
radius r™. 

We also did not follow the evolution of the bubble after it was cut off from the gap. 
Since it starts from a highly non-circular shape, it is expected to perform large amplitude 
oscillations. Remembering that the bubble is really a torus in three-dimensional space, 
it will also be unstable with respect to the Rayleigh instability ( Drazin fc Rcidj 1982) 
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Figure 5. The minimum radius rt as a function of time (dots), compared to the theoretical 
prediction ^J2~ro^/^\) (full line). Inset: the ratio r c /r c as a function of the initial radius varying 
between 1.25 • 10 -5 and 1.28 • 10~ 2 . The time for pinching r c was approximated using a linear 
extrapolation of z m in- The numerical results show very good agreement with the expected scaling 
law. 

and break up into a sequence of smaller bubbles. Evidently, this instability breaks the 
rotational symmetry and is thus well beyond the scope of the present work. 



5. Dynamics of retraction 

We now study the individual retraction events, characterized by a mass of fluid being 
accelerated by two line forces, in greater detail. Thus if 

dr b 

Itt = VUp 

is the velocity of the receding tip, the force balance reads 



- I M Up ^ 



dt 



dt 



(5.1) 



Inviscid coalescence of drops 



13 





t/r 3 b (6) ' " t/ r g(d) 

Figure 6. Two quantities characterizing retraction, Sr — rt(t) — r),(0) and z max as functions of 

time in rescaled units. Long-dashed and dotted lines represent power-law approximations to the 

early and long-time behavior, respectively. We find 8r oc t 2 for early times, while z max remains 



constant. For late times 8r oc t 0,8 and z ma x oc i ' 6 ; both behaviors are in agreement with (5.2). 
where M t i v is the mass being accelerated. This "added mass" is being pushed along by 
the structure of maximum radius z max that is forming at the end of the gap, and thus 



M t ip ~ Cz^ al . ( Landau fc Lifschitz 1982 ), where C is a numerical constant coming from 
the geometry of the void profile. Hence the equation of motion becomes 



dt \ max dt 



(5.2) 



For short times, the bubble does not have time to grow, so z max is approximately 
constant and given by the initial gap radius:z maa; « rj,(0)/2. This corresponds to a 
constant mass being accelerated by a constant force, and (|]^) leads to a quadratic 
growth of the retraction distance 5rb{t) = rb(t) — rb(0) oc t 2 . This is confirmed by the 
early time behavior of Sr^it) as shown in figure |6| Note that, consistent with ( pS.2[ ), z max 
remains constant. 

After this initial period of acceleration, the bubble radius z max starts to grow and the 
speed of retraction v t i P reaches a maximum, as seen in figure ^. This maximum must be 
set by the initial width w of the gap, and thus dimensional arguments lead to 



72A 



(5.3) 
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Figure 7. The speed of the retracting bridge vu p = drb(t)/dt as a function of time in rescaled 
units. The Culick- Taylor velocity (\/2 in these units) is represented by the dashed segment. 



The prefactor in (5.3) comes from balancing the inertial term v%/2 with the surface 
tension force k in (2.2), in analogy with the arguments of Culick and Taylor ( |Culick 



1960; Taylor 1959) for receding soap films. The curvature k has been approximated by 



1/w. As confirmed by figure ^, the maximum of Vu p is well approximated by the estimate 

After reaching a maximum, the speed of retraction decreases steadily, as the bubble 
grows and with it the added mass that has to be dragged along. The transversal bubble 
expansion is due to the rapid fluid motion along its sides which, according to Bernoulli's 
equation (|2.2| ), causes an under-pressure. Conversely, at the stagnation point behind the 
bubble the pressure is high and the bubble is curved inward (cf. figure [|). We do not yet 
have a fully quantitative theory of the bubble expansion, since this would require a precise 
knowledge of the bubble's shape. Namely, the fluid speed v m past the crest of the bubble 

VtipK c z max , in analogy to the flow 



is determined by its curvature k c (Lamb 1993): u 



past an ellipsoidal body. To close the system of equations, we would need an expression 
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for k c . However, we notice from figure ^| that the temporal growth of the bubble size 

Zmax is well described by a power law: z max oc t a e . Plugging this into equation (5.2) we 

find 

8r b <xt°- s , (5.4) 

in good agreement with simulations, cf. figure ^. The range of validity of the power laws 
proposed here can of course never exceed an order of magnitude, since the gap pinches 
off after time t ~ 10r^. 

Eventually, when the toroidal bubble separates from the gap, the velocity vu p has 
decreased to about half of v c . Therefore, the effect of the dynamical pressure v 2 ip /2 is 
reduced considerably relative to the capillary pressure. Numerically, we find that the 
capillary force is at least 4 times bigger than the dynamical pressure, which indicates 
that the velocity field can safely be neglected at reconnection, as we are forced to do 
owing to limitations of our numerical technique. 

6. Discussion 

We have shown that the merging of low viscosity fluid droplets leads to a self-similar 
sequence of void entrapments. It is interesting to note that the same power law behavior 
( |2.6| ) of rt, can be formally derived from a continuous evolution if vu p is assumed to be 



the Culick velocity (5.3). If the gap width w is estimated form the geometrical constraint 



r 2 , this immediately leads to c^n, m y/2/rb, which can be integrated to give a power 



t 1 / 2 . This is the argument given in (Eggers et al. 1999), which did not take reconnection 



into account. The reason it ends up to give the correct answer (apart from the prefactor) 
is that the size of the gap tip is rescaled to agree with the geometrical estimate (^]|) at 
each reconnection event. Thus although the bubble actually grows to a much larger size 
than r 2 , the balance implied by the above argument is actually true on average. 
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It might be equally tempting (Lister 2002| ) to apply the same reasoning to the force 

balance fl5.1|), by approximating (at least on average) the added mass by Mu p R* C^ nax oc 



rf. Integrating the corresponding equation of motion leads to r& oc t 2 ^ 5 . This apparent 
paradox is explained by the fact that the reconnection events destroy the momentum 
conservation implied by ( |5.l| ). Owing to bubble growth, momentum is distributed over a 
much larger volume than estimated from the simple geometrical argument. Accordingly, 
in the asymptotic limit of t <C 1 one obtains a motion that is faster than that given by 
the full calculation including reconnection. 

We would finally like to point out some questions inspired by this work. Firstly, it 
would be nice to develop a more complete theory of the bubble growth at the end of 
the receding meniscus. Secondly, we are not yet able to fully treat the velocity field 
after reconnection. Such a treatment may lead to an increase in fluctuations and perhaps 
some randomness during retraction. As pointed out in ( pguz fc Prospcrettj 1990), a finite 



velocity of approach will increase the likelihood of bubble entrapment during coalescence. 
Other interesting generalizations not yet considered in the present paper are the effect of 
an external fluid as well as viscous corrections. Clearly, a number of theoretical questions 
remain open. Perhaps more importantly, detailed experimental studies arc called for, for 
example to verify the phenomenon of bubble entrainment predicted by our analysis. 

It is our pleasure to thank Stephane Zaleski for its constant encouragement during this 
work. 
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